use "BP.formatted.dta", clear
set more off
local dv rep_civwar
local beta logGDP parliament disaster major_oil major_primary
local gamma logGDP parliament 

myoprob `dv' `beta'

preserve
drawnorm DH_b1-DH_b7, n(1000) means (e(b)) cov(e(V)) clear
save simulated_betas_op, replace
restore 
merge using simulated_betas_op
tab _merge
drop _merge
summarize DH*

postutil clear
postfile mfx p0 p1 p2 l0 u0 l1 u1 l2 u2 lgdpc using "lgdpmfxOP.dta", replace 
forvalues i = 5.8(.001)10.1 {

qui gen pr0 = normalden(-(DH_b1*`i'+DH_b6))*(-(DH_b1))
qui gen pr1 = (normalden(DH_b7-(DH_b1*`i'+DH_b6))-normalden(-(DH_b1*`i'+DH_b6)))*(-(DH_b1))
qui gen pr2 = normalden(DH_b7-(DH_b1*`i'+DH_b6))*(DH_b1)

sum pr0, meanonly
scalar _p0 = r(mean)
sum pr1, meanonly
scalar _p1 = r(mean)
sum pr2, meanonly
scalar _p2 = r(mean)

_pctile pr0, percentiles(2.5 97.5)
scalar _u0 = r(r2) 
scalar _l0 = r(r1)

_pctile pr1, percentiles(2.5 97.5)
scalar _u1 = r(r2) 
scalar _l1 = r(r1)

_pctile pr2, percentiles(2.5 97.5)
scalar _u2 = r(r2) 
scalar _l2 = r(r1)

post mfx (_p0) (_p1) (_p2) (_l0) (_u0) (_l1) (_u1) (_l2) (_u2) (`i')

drop pr0
drop pr1
drop pr2
}
postclose mfx

use "BP.formatted.dta", clear
set more off
local dv rep_civwar
local beta logGDP parliament disaster major_oil major_primary
local gamma logGDP parliament

ziop `dv' `beta', infl(`gamma')

preserve
drawnorm DH_b1-DH_b10, n(1000) means (e(b)) cov(e(V)) clear
save simulated_betas_z, replace
restore 
merge using simulated_betas_z
tab _merge
drop _merge
summarize DH*

postutil clear
postfile mfx p0 p1 p2 l0 u0 l1 u1 l2 u2 lgdpc using "lgdpmfx.dta", replace 
forvalues i = 5.8(.001)10.1 {
qui gen pr0 = (normal(-(DH_b6+`i'*DH_b1))-1)*normalden(DH_b9+`i'*DH_b7)*DH_b7-normal(DH_b9+`i'*DH_b7)*normalden(DH_b6+`i'*DH_b1)*DH_b1
qui gen pr1 = normal(-(DH_b6+`i'*DH_b1)-1)*normalden(DH_b9+`i'*DH_b7)*DH_b7-normal(DH_b9+`i'*DH_b7)*normalden(DH_b6+`i'*DH_b1)*DH_b1+( normalden(DH_b6+`i'*DH_b1)*normal(DH_b9+`i'*DH_b7)-normalden(-(DH_b6+`i'*DH_b1))*normal(DH_b9+`i'*DH_b7) )*DH_b1
qui gen pr2 = (normal((DH_b10-(DH_b6+`i'*DH_b1)))-normal(-(DH_b6+`i'*DH_b1)))*normalden(DH_b9+`i'*DH_b7)*DH_b7+( normalden(-(DH_b6+`i'*DH_b1))*normal(DH_b9+`i'*DH_b7)-normalden(DH_b10-(DH_b6+`i'*DH_b1))*normal(DH_b9+`i'*DH_b7) )*DH_b1

sum pr0, meanonly
scalar _pr0 = r(mean)
sum pr1, meanonly
scalar _pr1 = r(mean)
sum pr2, meanonly
scalar _pr2 = r(mean)

_pctile pr0, percentiles(2.5 97.5)
scalar _u0 = r(r2)
scalar _l0 = r(r1)
_pctile pr1, percentiles(2.5 97.5)
scalar _u1 = r(r2)
scalar _l1 = r(r1)
_pctile pr2, percentiles(2.5 97.5)
scalar _u2 = r(r2)
scalar _l2 = r(r1)

post mfx (_pr0) (_pr1) (_pr2) (_l0) (_u0) (_l1) (_u1) (_l2) (_u2) (`i')
drop pr0
drop pr1
drop pr2
}
postclose mfx

use "BP.formatted.dta", clear
set more off
local dv rep_civwar
local beta logGDP parliament disaster major_oil major_primary
local gamma logGDP parliament

ziopc `dv' `beta', infl(`gamma') 

preserve
drawnorm DH_b1-DH_b11, n(1000) means (e(b)) cov(e(V)) clear
save simulated_betas_zc, replace
restore 
merge using simulated_betas_zc
tab _merge
drop _merge
summarize DH*

postutil clear
postfile mfx p0 p1 p2 l0 u0 l1 u1 l2 u2 lgdpc using "lgdpmfxC.dta", replace 
forvalues i = 5.8(.001)10.1 {
qui gen pr0 = (normal((-(DH_b6+`i'*DH_b1)+(DH_b11*(DH_b9+`i'*DH_b7)))/sqrt(1-(DH_b11*DH_b11)))-1)*normalden(DH_b9+`i'*DH_b7)*DH_b7-normal(((DH_b9+`i'*DH_b7)-(DH_b11*(DH_b6+`i'*DH_b1)))/sqrt(1-(DH_b11*DH_b11)))*normalden(DH_b6+`i'*DH_b1)*DH_b1
qui gen pr1 = (normal((DH_b10-(DH_b6+`i'*DH_b1)+(DH_b11*(DH_b9+`i'*DH_b7)))/sqrt(1-(DH_b11*DH_b11)))-normal((-(DH_b6+`i'*DH_b1)+(DH_b11*(DH_b9+`i'*DH_b7)))/sqrt(1-(DH_b11*DH_b11))))*normalden(DH_b9+`i'*DH_b7)*DH_b7+( normalden(-(DH_b6+`i'*DH_b1))*normal(((DH_b9+`i'*DH_b7)+(DH_b11*(-(DH_b6+`i'*DH_b1))))/sqrt(1-(DH_b11*DH_b11)))-normalden(DH_b10-(DH_b6+`i'*DH_b1))*normal(((DH_b9+`i'*DH_b7)+(DH_b11*(DH_b10-(DH_b6+`i'*DH_b1))))/sqrt(1-(DH_b11*DH_b11))))*DH_b1
qui gen pr2 = normal(((DH_b6+`i'*DH_b1)-(DH_b10)-(DH_b11*(DH_b9+`i'*DH_b7)))/sqrt(1-(DH_b11*DH_b11)))*normalden(DH_b9+`i'*DH_b7)*DH_b7+normalden((DH_b6+`i'*DH_b1)-DH_b10)*normal(((DH_b9+`i'*DH_b7)-(DH_b11*((DH_b6+`i'*DH_b1)-DH_b10)))/sqrt(1-(DH_b11*DH_b11)))*DH_b1

sum pr0, meanonly
scalar _pr0 = r(mean)
sum pr1, meanonly
scalar _pr1 = r(mean)
sum pr2, meanonly
scalar _pr2 = r(mean)

_pctile pr0, percentiles(2.5 97.5)
scalar _u0 = r(r2)
scalar _l0 = r(r1)
_pctile pr1, percentiles(2.5 97.5)
scalar _u1 = r(r2)
scalar _l1 = r(r1)
_pctile pr2, percentiles(2.5 97.5)
scalar _u2 = r(r2)
scalar _l2 = r(r1)

post mfx (_pr0) (_pr1) (_pr2) (_l0) (_u0) (_l1) (_u1) (_l2) (_u2) (`i')
drop pr0
drop pr1
drop pr2
}
postclose mfx

cd "C:\Users\Danny\Documents\ZIOP\ZiOP_Datasets\ZiOP Datasets"
use "lgdpmfxOP.dta", clear

twoway (line p0 lgdpc, lpattern(solid)) (line l0 lgdpc, lpattern(dash)) (line u0 lgdpc, lpattern(dash)), yline(0) xtitle(ln(GDP/capita)) legend(off) title(OP - Pr(No Violence)) scheme(s2mono)
graph save 0OP, replace

twoway (line p1 lgdpc, lpattern(solid)) (line l1 lgdpc, lpattern(dash)) (line u1 lgdpc, lpattern(dash)), yline(0) xtitle(ln(GDP/capita)) legend(off) title(OP - Pr(Repression)) scheme(s2mono)
graph save 1OP, replace

twoway (line p2 lgdpc, lpattern(solid)) (line l2 lgdpc, lpattern(dash)) (line u2 lgdpc, lpattern(dash)), yline(0) xtitle(ln(GDP/capita)) legend(off) title(OP - Pr(Civil War)) scheme(s2mono)
graph save 2OP, replace

graph combine 0OP.gph 1OP.gph 2OP.gph, ycommon r(1)

graph save OP, replace

graph export lgdpOP.eps, replace
graph export lgdpOP.png, replace

use lgdpmfx.dta, clear

twoway (line p0 lgdpc, lpattern(solid)) (line l0 lgdpc, lpattern(dash)) (line u0 lgdpc, lpattern(dash)), yline(0) xtitle(GDP/capita) legend(off) title(ZiOP - Pr(No Violence)) scheme(s2mono)
graph save 0Z, replace

twoway (line p1 lgdpc, lpattern(solid)) (line l1 lgdpc, lpattern(dash)) (line u1 lgdpc, lpattern(dash)), yline(0) xtitle(GDP/capita) legend(off) title(ZiOP - Pr(Repression)) scheme(s2mono)
graph save 1Z, replace

twoway (line p2 lgdpc, lpattern(solid)) (line l2 lgdpc, lpattern(dash)) (line u2 lgdpc, lpattern(dash)), yline(0) xtitle(GDP/capita) legend(off) title(ZiOP - Pr(Civil War)) scheme(s2mono)
graph save 2Z, replace

graph combine 0Z.gph 1Z.gph 2Z.gph, ycommon r(1)

graph save Z, replace

graph export lgdpZ.eps, replace
graph export lgdpZ.png, replace

use lgdpmfxC.dta, clear

twoway (line p0 lgdpc, lpattern(solid)) (line l0 lgdpc, lpattern(dash)) (line u0 lgdpc, lpattern(dash)), yline(0, lcolor(black)) xtitle(ln(GDP/capita)) legend(off) title(ZiOPC - Pr(No Violence)) scheme(s2mono)
graph save 0ZC, replace

twoway (line p1 lgdpc, lpattern(solid)) (line l1 lgdpc, lpattern(dash)) (line u1 lgdpc, lpattern(dash)), yline(0) xtitle(ln(GDP/capita)) legend(off) title(ZiOPC - Pr(Repression)) scheme(s2mono)
graph save 1ZC, replace

twoway (line p2 lgdpc, lpattern(solid)) (line l2 lgdpc, lpattern(dash)) (line u2 lgdpc, lpattern(dash)), yline(0) xtitle(ln(GDP/capita)) legend(off) title(ZiOPC - Pr(Civil War)) scheme(s2mono)
graph save 2ZC, replace

graph combine 0ZC.gph 1ZC.gph 2ZC.gph, ycommon r(1)

graph save ZC, replace

graph export lgdpC.eps, replace
graph export lgdpC.png, replace

graph combine OP.gph Z.gph ZC.gph, ycommon r(3)

graph export lgdpAll.eps, replace
graph export lgdpAll.png, replace